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ABSTRACT 

Determining the virial factor of the broad-line region (BLR) gas is crucial for calibrating AGN black hole 
mass estimators, since the measured line-of-sight velocity needs to be converted into the intrinsic virial ve- 
locity. The average virial factor has been empirically calibrated based on the M B H-er» relation of quiescent 
galaxies, but the claimed values differ by a factor of two in recent studies. We investigate the origin of the 
difference by measuring the M BH - <t* relation using an updated galaxy sample from the literature, and explore 
the dependence of the virial factor on various fitting methods. We find that the discrepancy is primarily caused 
by the sample selection, while the difference stemming from the various regression methods is marginal. How- 
ever, we generally prefer the FITEXY and Bayesian estimators based on Monte Carlo simulations for the 
Msh-(J* relation. In addition, the choice of independent variable in the regression leads to ~ 0.2 dex variation 
in the virial factor inferred from the calibration process. Based on the determined virial factor, we present the 
updated M BH - c* relation of local active galaxies. 

Subject headings: galaxies: quiescent - galaxies: active - scaling relation - virial factor - regression methods 



1. INTRODUCTION 

Supermassive black holes (SMBHs) are thought to be ubiq- 
uitous in the centers of virtually all massive galaxies (e.g., Ko- 
rmendy & Richstone 1995; Richstone et al. 1998; Ferrarese 
& Ford 2005). The close connection of black hole growth to 
galaxy evolution is inferred from the discovery of tight corre- 
lations between the masses of SMBHs (M B h) and the global 
properties of host galaxies, such as the stellar velocity disper- 
sion (cr»; Ferraresse & Merritt 2000; Gebhardt et al. 2000) 
bulge luminosity (Lbui; Magorrian et al. 1998; Marconi & 
Hunt 2003), and bulge mass (M bul ; Haring & Rix 2004). The 
origin of these connections has been investigated in theoret- 
ical studies of galaxy evolution either through the introduc- 
tion of active galactic nucleus (AGN) feedback (e.g., Kauff- 
mann & Haehnelt 2000; Di Matteo et al. 2005; Croton et al. 
2006; Hopkins et al. 2006; Bower et al. 2006; Somerville 
et al. 2008; Booth & Schaye 2009) or as simply being the 
result of a hierarchical merging framework (e.g., Peng 2007; 
Hirschmann et al. 2010; Janke & Maccio 2011). The in- 
terplay between the black holes and galaxies is now one of 
the basic ingredients in our understanding of galaxy forma- 
tion and evolution. 

In order to better understand the origin and evolution of 
the SMBH-host galaxy connection, AGN demographics, and 
the growth of the SMBHs through cosmic time, an accu- 
rate and precise measurement of black hole mass is essen- 
tial. Stellar/gas dynamical modeling is commonly used to 
measure the black hole masses in quiescent galaxies. How- 
ever, this technique requires high spatial resolution to re- 
solve the sphere-of-influence of the black hole, thereby lim- 
iting it to the local universe. In active galaxies, the black 
hole mass can be determined by utilizing AGN variability. 
The reverberation mapping technique (Blandford & McKee 
1982; Peterson 1993) has been used to estimate the mean 
size of the broad line region (BLR, Rblr) by cross-correlating 
the continuum light curve with the broad emission line light 
curve. Combining Rblr with the line-of-sight velocity width 
(AV) measured from the variable component of the broad 
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emission line provides a virial black hole mass estimate as 
M B h = f AV 2 Rblr/G, where G is the gravitational constant 
and / is the virial factor that converts the measured virial 
product into the actual black hole mass. This technique is 
also limited to around 50 AGNs to date since it requires exten- 
sive photometric and spectroscopic monitoring observations. 
This technique has established the empirical size-luminosity 
relation (Wandel et al. 1999; Kaspi et al. 2000; Bentz et 
al. 2006a, 2009a), which is the basis for the single-epoch 
(SE) method. In the SE method, one simply substitutes the 
time-consuming BLR size measurement with AGN luminos- 
ity using the size-luminosity relation. This therefore provides 
estimates of black hole masses for broad line AGNs from a 
single spectroscopic observation, thus expanding the sample 
size substantially at any redshift. However, both methods suf- 
fer from the large uncertainty stemming from the unknown 
virial factor (see Park et al. 2012), which depends on the ge- 
ometry and kinematics of BLR of individual AGNs. 

Instead, an empirically calibrated average virial factor has 
been applied to most AGN black hole mass estimates, ex- 
cept for only a few objects where dynamical mass measure- 
ments can be obtained (e.g., Davies et al. 2006; Onken et 
al. 2007; Hicks & Malkan 2008). The first calibration of the 
virial factor was performed by Onken et al. (2004). They 
derived (/) = 5.5 ±1.8 based on a sample of 14 AGNs, for 
which both reverberation masses and stellar velocity disper- 
sions were available, by forcing the AGN host galaxies to 
obey the same Mbh-u* relationship as for quiescent galaxies. 
By enlarging the dynamical range of the AGN sample, Woo 
et al. (2010) determined the virial factor as log(/) = 0.72+g % 
(i.e., (/) = 5.2 ± 1 .2) based on an updated reverberation sam- 
ple of 24 AGNs, which included 8 low-mass local Seyfert 1 
galaxies from the Lick AGN Monitoring project (Bentz et al. 
2009b). They provided the upper limit of uncertainty in the 
derived virial factor as 0.43 dex based on the intrinsic scat- 
ter in the relation. In contrast, Graham et al. (201 1) reported 
(/) = 2.8^} 5, based on their updated M B h-<t» relation of qui- 
escent galaxies, and an upated AGN sample, which is a factor 
of 2 smaller than the aforementioned values. Graham et al. 
(201 1) commented that the value of (/) might be even further 
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lowered due to the effect of radiation pressure (see Marconi 
et al. 2008). This correspondingly reduces the black hole 
mass estimates for most AGNs by that amount, influencing 
all of the studies incorporating single-epoch AGN black hole 
masses. Thus it is important to investigate the origin of this 
difference and check for possible biases in the calibration pro- 
cess. 

Since the derived Mbh - relation of quiescent galaxies is 
used to calibrate the virial factor in AGN mass estimators, un- 
der the assumption that the same Mbh _ a* relation holds for 
AGN host galaxies, it is important to investigate the differ- 
ences in the Mbh _ c* relations of quiescent galaxies reported 
in the literature, and to study their effect on the derived virial 
factors. Originally the slopes of theMBH-cr* relation reported 
by Ferrarese & Merritt (2000) and Gebhardt et al. (2000) were 
4.8 ±0.5 and 3.75 ±0.3 based on 12 and 26 galaxy samples, 
respectively. After that, various slopes have been reported 
in the literature, ranging from 3.68 to 5.95. Although the 
slopes are roughly consistent with the theoretical expectations 
of M oc cr 5 (Silk & Rees 1998) and M oc a 4 (Fabian 1999), 
their difference and change are noteworthy. The possible ori- 
gin of the difference in slopes has been investigated in the 
literature. The related factors are: (1) the type of regression 
method adopted (Tremaine et al. 2002; Novak et al. 2006; see 
Kelly 2007 for general applications of regression), (2) the size 
of the assigned uncertainty on the velocity dispersion (Merritt 
& Ferrarese 2001; Tremaine et al. 2002), (3) the velocity dis- 
persion measures used (Tremaine et al. 2002), (4) the adopted 
value of velocity dispersion for the Milky Way (Merritt & Fer- 
rarese 2001; Tremaine et al. 2002), (5) the spatial resolution 
of the data for the resolved BH sphere-of-influence (Ferrarese 
& Ford 2005; see also Gultekin et al. 2009, 2011; Batchel- 
dor 2010) (6) the morphological type of the sample used (Hu 
2008; Graham 2008; Gultekin et al. 2009; Greene et al. 2010; 
Graham et al. 2011). 

To understand the origin of the differences in the derived 
Msn-cr* relationships, we investigate in this work 3 main 
issues: the difference in samples, the difference in regres- 
sion methods, and the direction of the regression analysis. 
Recently, with new measurements and improved modeling 
the number of dynamical mass measurements is continuously 
growing both at the high-mass and low-mass end regimes. To 
date, a total of 67 black hole masses in quiescent galaxies has 
been measured via stellar/gas/maser kinematics (see the most 
recent compilation from McConnell et al. 2011 and refer- 
ences therein). Therefore, it is presently a good time to inves- 
tigate what effect the difference in samples has on the derived 
Mbh _ c» relation using the largest sample ever. In addition, 
various estimators have been used for the regression analy- 
sis in the black hole scaling relation studies: FITEXY (e.g., 
Tremaine et al. 2002, Novak et al. 2006, Kim et al. 2008, Li 
et al. 2011, Beifiori et al. 2012, Vika et al. 2012, McConnell 
et al. 2011), BCES (e.g., Ferrarese & Merritt 2000, Ferrarese 
& Ford 2005, Hu 2008, Bentz et al. 2009a, Bennert et al. 

2010, Graham et al. 201 1), Maximum likelihood (e.g., 
Gultekin et al. 2009, Greene et al. 2010, Schulze & Gebhardt 
2011), Bayesian approach (linmix_err) (e.g., Sani et al. 

201 1, Xiao et al. 201 1, Mancini & Feoli 2012). Thus, in or- 
der to investigate the difference in the derived scaling relation- 
ships caused by the sample selection, it is important to inves- 
tigate differences between the estimators for the Mbh _ c* re- 
lation analysis. Finally, adopting the choice of independent 
variable is another issue for determining the Mbh _ c* rela- 
tion. Motivated by the suggestion by Graham et al. (201 1) to 
use the 'inverse' fit to calibrate the single-epoch AGN mass 
estimates, we present results based on both of the forward and 



inverse regressions. 

This paper is organized as follows. In the next section, we 
describe the most commonly used regression methods in as- 
tronomy with their explicit implementations. In Section [3] 
we re-measure the M B h - relation using 3 different sam- 
ples from the literature and investigate the difference due to 
the regression methods and samples. In Section|4]we present 
our main result for the calibrated virial factors and discuss 
the difference based on the regression methods and samples. 
The difference from the direction of regression is discussed in 
Section|5] Finally, we summarize and conclude in Section|6] 

2. LINEAR REGRESSION TECHNIQUES 

Linear regression methods^in astronomy were exhaustively 
discussed in the pioneering paper, Isobe et al. (1990). They 
provided formulae for 5 unweighted bivariate linear regres- 
sion coefficients with their error estimates, and recommended 
the bisector line for the case of treating the variables sym- 
metrically. The second paper in the series, Feigelson & 
Babu (1992), extended their work by accommodating boot- 
strap and jackknife resampling procedures for error estima- 
tion, weighted regression, and truncated/censored regressions. 
In addition, they suggested practical strategies for linear re- 
gression problems in astronomy. 

The BCES estimator (Bivariate Correlated Errors and in- 
trinsic Scatter) was proposed by Akritas & Bershady (1996) 
in order to incorporate, heteroscedastic measurement errors, 
intrinsic scatter, and correlation in the measurement errors. 
The method of minimizing a x 2 statistic (FITEXY), which 
account for measurement error in both the dependent and in- 
dependent variable, was modified by Tremaine et al. (2002) 
to incorporate intrinsic scatter. They added the unknown con- 
stant intrinsic variance term in quadrature to the error of the 
dependent variable and determined it so that the reduced x 2 
is equal to a value of unity. Based on the Monte Carlo simu- 
lations performed by Tremaine et al. (2002) and Novak et al. 
(2006), they concluded that the modified FITEXY is a better 
estimator than the BCES. In particular, Tremaine et al. (2002) 
concluded that the BCES tends to be biased when the sample 
size is small or the mean square of the x errors is comparable 
to the variance of x distribution, and that it becomes ineffi- 
cient when there is a single measurement with much larger 
error than others. 

Kelly (2007) developed a sophisticated Bayesian linear re- 
gression technique, termed linmix_err. It accounts for 
intrinsic scatter in the relationship, heteroscedastic measure- 
ment errors in both the independent and dependent variables, 
and correlation between the measurement errors. This method 
uses a Gaussian mixture model for the distribution of inde- 
pendent variables, which is shown to work well particularly 
when the measurement errors are large by avoiding the bias 
incorporated if the choice of x-distribution model is incorrect 
(also noted in Auger et al. 2010). The method assumes that 
the measurement errors and intrinsic scatter are Gaussian, and 
it accommodates multiple independent variables, nondetec- 
tions, and selection effects. 

Recently, Gultekin et al. (2009) applied a maximum 
likelihood method to determine the M - a and M-L rela- 
tions by naturally incorporating an intrinsic scatter and upper 
limits. They also extensively investigated the distributional 
forms for the measurement error and intrinsic scatter. How- 
ever, they did not include a model for the distribution of the 
independent variable, but instead used Monte Carlo sampling 
to assess the impact of the measurement errors in the indepen- 

2 For recent reviews, please see Hogg et al. (2010) and Caimmi (201 la,b). 
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dent variable on the parameter estimates. 

To sum up, the four methods for linear regression that 
have been used to characterize the black hole/host galaxy 
scaling relationships are: (1) BCES (Akritas & Bershady 
1996), (2) FITEXY (Tremaine et al. 2002), (3) Maximum 
likelihood (Giiltekin et al. 2009), (4) LINMIX_ERR 
(Kelly 2007). In this section we explicitly show our imple- 
mentation and usage of each method. We assume the model 
of y = a + fix in the following analysis. 

2.1. BCES estimator 

The BCES(F|X) estimator is implemented using the for- 
mula described in Akritas & Bershady (1996), 



P 



cov(x,y)-(a xy ) 



var(x) - {a]) 

a =(y)-P( x ), 



(1) 
(2) 



where co\(x,y) is the covariance of x and y, a x is the standard 
deviation of the measurement error (i.e., standard measure- 
ment error) in x, var(x) is the variance of x, and a xy is the 
covariance between the measurement errors in x and y. The 
intrinsic variance (i.e., variance in the intrinsic scatter) is esti- 
mated following the expression given in Cheng & Riu (2006) 
and Kelly (2007), 



:= \J var(y) - (of - f3 [cov(x,y) - {a xy )\ . 



(3) 



The uncertainties in the parameters can be estimated with the 
bootstrap or using analytical estimates given in Akritas & Ber- 
shady (1996). In this work we assume a xy = (i.e., uncorre- 
cted measurement errors), as most values of x and y in the 
A/jjh - cr* samples were independently measured and the co- 
variances between the measurement errors are not provided 
in the literature. Thus, simply assuming the zero covariance 
seems to be more reasonable for these very heterogeneously 
collected Mbh - c* samples. In addition, there is no result 
incorporating the correlated measurement errors (i.e., a^) to 
M%K-a* fitting in the literature at least to our knowledge. 
Thus, to compare consistently with the results from the liter- 
ature we here set a xy = 0. 

2.2. FITEXY estimator 

The FITEXY (Press et al. 1992), modified by Tremaine 
et al. (2002), is implemented in our work in IDL using 
the mpfit (Markwardt 2009) Levenberg-Marquardt least- 
squares minimization routine. Note that our implementation 
is basically similar to that given in Williams et al. (2010)0. It 
performs the linear regression by minimizing 
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where a and /3 are the regression coefficients, a x and a y are 
the standard deviation in the measurement errors, and af nt is 
the intrinsic variance. The value of a mi is iteratively adjusted 
as an effective additional y error by repeating the fit until one 
obtains x 2 /(N-2) = 1 (i.e., following the suggested iterative 
procedure given in Bedregal et al. 2006 and Bamford et al. 
2006). If after the initial iteration the reduced \ 2 is l ess than 
one, then no further iterations occur and one sets a m = 0. We 
estimate uncertainties in the regression parameters with the 
bootstrap method. 



2.3. Maximum Likelihood estimator 

The method of maximum likelihood is implemented 
similarly as given in Giiltekin et al. (2009) (see also Woo et 
al. 2010). Under the assumptions of uncorrelated Gaussian 
measurement errors in both coordinates and Gaussian intrin- 
sic scatter along y, the Gaussian likelihood function is given 

by 



i 



exp 



(y,— a-/3xj) 
2a? 



where 



2 2 
cr,- =0- y<i + l 



j2 2 ,2 
i cr,., + (T int 



(5) 



(6) 



Note that this likelihood function implicitly assumes that 
the independent variable is uniformly distributed (i.e., a uni- 
form prior for the intrinsic distribution of x). Then the log- 
likelihood function is 



-21n£ = ^ln (27rer, 2 )+^] 

N 



2 



(7) 



This likelihood approach is more complete than the \ 2 
method given in Equation (J4]i in a sense that it includes the in- 
trinsic variance term in both the normalization and exponent 
of the likelihood function, and determines it simultaneously 
with the other regression coefficients. To estimate the best-fit 
parameters of (a,f3,ai nt ) using the maximum likelihood es- 
timate (MLE), we minimize -21n£ using the downhill sim- 
plex method implemented as AMOEBA (Press et al. 1992) in 
IDL. Then we adopt uncertainties in parameters as the average 
difference where In C decreases from its maximum value by 
0.5. We also estimate the parameter errors using the bootstrap 
method. 



2.4. Bayesian estimator (Unmix _err) 

The Bayesian linear regression routine, linmix_err, de- 
veloped by Kelly (2007) is available from the NASA IDL 
astronomy user's librarjy. Here we briefly summarize the 
method. For details, please refer to Kelly (2007). 

This method assumes Gaussian intrinsic scatter, Gaussian 
measurement errors, and a weighted sum of K Gaussian func- 
tions for the distribution of the independent variable. The 
choice of a Gaussian mixture model was motivated in that it 
can not only approximate well various intrinsic distributions 
of the independent variable, but it is also a mathematically 
convenient conjugate family. The full measured data likeli- 
hood function is expressed as a mixture of bivariate normal 
distributions 



ll\ 



i=l k=l 



~izi-Ck) T V k -j(zi-Ck) 



, (8) 



where J2k=\ = 1 an d the measured data, means, and covari- 



3 http://purl.org/mike/mpfitexy 
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ance matrices are, respectively, 



-f2 



Ck-- 



V k4 = 



/3r A 2 + cr.v V ,,- Tl + a 2 u 



(9) 
(10) 
(11) 



In order to calculate the posterior probability distribution 
of the model parameters for the given measured data, it 
adopts uniform prior distributions on the regression param- 
eters (a,P,af nt ). It also adopt a Dirichlet, normal, and scaled 
inverse-^ 2 prior on the mixture model parameters (n^, /ik, t\), 
respectively. The data is 'fit' using a Markov Chain Monte 
Carlo (MCMC) sampler. For each regression parameter, we 
take the best-fit value and uncertainty as the posterior me- 
dian and posterior standard deviation from the marginal pos- 
terior distributions using the 100,000 random draws returned 
by the MCMC sampler. Note that the likelihood function 
given in Equation (0 converges to that of the maximum 
likelihood method given in Equation (0 if the distribu- 
tion for the independent variable is assumed to be uniform 
rather than a mixture of normals and the measurement errors 
are uncorrected. 

As an illustration, we also used the likelihood function in 
the case of a single Gaussian model (K = l,%k = 1) for the 
distribution of independent variable with uncorrected mea- 
surement errors (<r xy ,i = 0). For comparison with the procedure 
assu min g a uniform intrinsic distribution (MLE) given in Sec- 
tion [53] we compute the maximum likelihood estimate (i.e., 
MLEig) utilizing the likelihood function derived from assum- 
ing the distribution of independent variable is a Gaussian. 

In the following sections, we determine the a, (3, and a mt 
parameters with the corresponding error estimates for the 
Mbh - cr* relation using each regression technique described 
above. 

3. THE M BH -<?* RELATIONS 

The Mbh _ c* relation is generally expressed as the log- 
linear form, 

log(M BH /M ) = a + /31og(cr»/200 km s" 1 ). (12) 

Here y = log(M B H/M Q ) and x = log(cr»/200 kms" 1 ). 
We assume that the measurement errors in the log- 
arithms of mass and stellar velocity dispersion are 
symmetric by taking the symmetric interval in log 
space, i.e., ei og M BH = (log M^ 61 -log M^™ 61 ") /2 and ei ogo -, = 
(loger" ppel -logo-J° wer ) /2. Following Graham et al., for 
their data we assume that the measurement errors on 
the logarithm of mass are symmetric by taking the av- 
erage of upper and lower la uncertainties on the lin- 
ear scale and propagating it onto the logarithmic scale, 



ClogMeH 



= 0.5 (( 



MbT + ^hO/^bh^IO). The measure- 



ment errors in the logarithm of the dispersion are ei ogCTi> = 



0.5 (< 



upper + ^lower^ 



/ (er» In 10). However, note that this choice 



of error bars does not significantly affect our results. 

3.1. Re-measuring the relation with four methods 

In this section we consistently re-derive the Mbh _ c* rela- 
tion using three literature samples to check our implementa- 
tion of the fitting methods. Figure [T] shows the re-estimated 
Mbh-c* relations of the data from Gultekin et al. 2009 (top), 
McConnell et al. 2011 (middle), and Graham et al. 2011 



O 



o 



m 8 



o 



OLS 

FITEXY 
BCES 
Bayesian 
MLE 




log (<J„/km s ') 



FIG. 1. — The Mbh - relations using each dataset from Gultekin et 
al. 2009 (top), McConnell et al. 2011 (middle), and Graham et al. 2011 
(bottom). Each regression line is derived from five different methods (see 
the text and Table 0. The sample distributions for the logarithms of black 
hole masses and stellar velocity dispersions are shown in the right side and 
top side of each panel as grey histograms. The non-overlapping sample in 
between the McConnell and Graham data are marked with a filled dot inside 
open circles. Note that the common sample has 50 galaxies. The barred 
galaxies are marked with an open square enclosing open circles. 
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TABLE 1 

Re-estimation of Parameters for the M B h -f, Relation of Quiescent Galaxy Samples: 
log(M BH /Mq ) = a + log(<7» /200 km s" 1 ) 



Method 




Forward regression 






Inverse regression" 








P 




a = -a im /(S im 


/3=l/ftnv 


°"int — °"int.inv//^inv 








Gultekin et al. (2009) Sample b 






OLS 


8.18 ±0.06 


3.91 ±0.28 




8.21 ±0.07 


5.60 ±0.68 




MLE 


8.19±0.06 


3.97 ±0.31 


0.39 ±0.06 


8.22 ±0.07 


5.61 ±0.70 


0.47 ±0.08 


BCES 


8.18±0.06 


4.01 ±0.32 


0.38 ±0.07 


8.20 ±0.07 


5.27 ±0.70 


0.43 ±0.10 


FITEXY 


8.19±0.06 


4.06 ±0.32 


0.39 ±0.06 


8.21 ±0.07 


5.35 ±0.66 


0.45 ±0.09 


Bayesian 


8.19±0.07 


4.04 ±0.40 


0.42 ± 0.05 


8.21 ±0.08 


5.44 ±0.56 


0.49 ±0.09 








McConnell 


etal. (2011) Sample 1 - 







OLS 


8.27 ±0.06 


4.87 ±0.37 




8.33 ±0.06 


6.55 ±0.50 




MLE 


8.28 ±0.06 


4.92 ±0.34 


0.41 ±0.05 


8.33 ±0.06 


6.43 ±0.51 


0.47 ±0.06 


BCES 


8.28 ±0.06 


5.06 ±0.41 


0.43 ± 0.05 


8.33 ±0.06 


6.36±0.51 


0.48 ±0.07 


FITEXY 


8.28 ±0.06 


5.07 ±0.36 


0.43 ± 0.05 


8.32 ±0.06 


6.29 ±0.49 


0.47 ±0.06 


Bayesian 


8.27 ±0.06 


5.06 ±0.36 


0.44 ± 0.05 


8.32 ±0.07 


6.31 ±0.46 


0.49 ±0.07 



Graham et al. (2011) Sample 



OLS 

MLE 

BCES 

FITEXY 

Bayesian 



.13±0.05 
.14 ±0.05 
.13±0.05 
.15±0.05 
.15±0.05 



4.75 ±0.29 
4.72 ±0.29 
5.13 ±0.35 
5.08 ±0.34 
5.08 ±0.36 



0.30 ±0.04 
0.31 ±0.04 
0.31 ±0.04 
0.31 ±0.05 



3. 16 ±0.06 

3.17 ±0.06 
3.15 ±0.06 
3.16±0.05 
3. 17 ±0.06 



6.22 ±0.46 
6.06 ±0.46 
5.95 ±0.45 

5.84 ±0.42 

5.85 ±0.42 



0.33 ±0.05 
0.34 ±0.05 
0.33 ±0.05 
0.34 ±0.06 



NOTE. — Forward regression=fit logMgH on logo - *; Inverse regression=fit logo - * on logMBH; OLS=Ordinary Least Squares; 
MLE=Maximum Likelihood Estimates; BCES=estimator of Akritas & Bershady (1996); FITEXY=estimator of Tremaine et al. 
(2002); Bayesian=Bayesian posterior median estimates using linmix_err procedure of Kelly (2007). 
a Inverse regression and its results will be discussed in Section|3] 

b We used 49 galaxies listed in Table 1 in Gultekin et al. (2009) without upper limits for comparison with other samples. For the 
MLE estimate, we also estimate the parameters using the same likelihood function and error estimation method given in Gultekin et 
al. (2009). The result is («, 0, cr int ) = (8. 18 ± 0.06, 3.97 ± 0.39, 0.42 ± 0.05). Note that there are two different mass measurements for 
NGC 1399 and NGC 5128. 

c We used all of the 67 galaxies listed in Table 4 in McConnell et al. (201 1), while they used only 65. We found that there is a typo in 
the M B h of NGC 1023 in their Table 4. Thus we corrected the value from 14.6 X 10 7 to 4.6 X 10 7 . Note that there are two different 
mass measurements for the NGC 1399 and for the NGC 5128. Following their scheme, if we apply half weights for them, then we 
get the same result with that of McConnell etal. (2011), i.e., («, /3,cr mt ) = (8.28 ± 0.06,5.13 ± 0.34,0.42 ± 0.05) for the FITEXY 
estimate. 



(bottom) using the various methods described in the previous 
section. We also include the ordinary least squares (hereafter 
OLS) line as a reference. For the sample of 49 galaxies with- 
out upper limits from Table 1 in Gultekin et al. (2009), we 
also follow the same fitting scheme of their maximum likeli- 
hood estimator, which is slightly different to the method de- 
scribed in Section|23] Using the Gultekin et al. (2009) proce- 
dure, we first perform the fit without accounting for the mea- 
surement errors in the independent variables. Then, we incor- 
porate the effects of measurement errors in x into the param- 
eter uncertainties by adding in quadrature the standard devia- 
tions estimated from the Monte Carlo fitting results for 10 4 re- 
alizations of the independent variables. Recently, McConnell 
et al. (201 1) have updated the compiled sample of 49 galaxies 
from Gultekin et al. (2009) by including new black hole mass 
measurements and revising earlier black hole masses based 
on improved stellar orbit modeling, which accounts for dark 
matter halos (Gebhardt & Thomas 2009; van den Bosch & 
de Zeeuw 2010; Shen & Gebhardt 2010; Schulze & Geb- 
hardt 2011). We have consistently re-estimated slopes of the 
Mbh - C* relation using the 67 galaxies listed in Table 4 of 
McConnell et al. (2011). The independently compiled sam- 
ple of 64 galaxies from Graham et al. (2011) is also used. 
TableQ]lists all regression results. Note that we get consistent 
results with each paper if we choose the same method and 
setting used by the respective papers. 

For the Gultekin et al. sample the difference between fit- 
ted lines is only marginal since they assumed a minimum of 
5% measurement errors on er*; such small errors in a* are 



found to have relatively small impact on the regression result, 
as described in the next section. However, the slope of the 
Mbh - &* relation derived from the updated sample of Mc- 
Connell et al. is significantly larger than that of the Gultekin 
et al., increasing from ~ 4 to ~ 5. As implied by the his- 
tograms for the black hole mass distributions shown in grey 
in the figures, this increase is mostly due to inclusion of both 
the low-mass and high-mass sample, which generally show an 
offset trend to the relation of the Gultekin et al. sample. For 
the Graham et al. sample the difference between fitted lines is 
marginally significant, and the slope is apparently divided into 
two groups, since they assigned relatively large (10%) mea- 
surement errors in er*. It seems that the bias of the maximum 
likelihood estimator starts to become non-negligible for 
x-errors of this magnitude. Therefore we investigate in details 
the effect of the amplitude of the x-error on these estimators 
in the following sections. 

3.2. The effect of the adopted measurement uncertainty of 

Merritt & Ferrarese (2001) first noted that ignoring the mea- 
surement errors in the velocity dispersion leads to a biased 
slope (i.e., underestimates) in the Mbh _ c* relation. How- 
ever, Tremaine et al. (2002) argued that the effect of the mea- 
surement errors in the velocity dispersion is not significant 
even at the 10% error level using two estimators, BCES and 
FITEXY, based on their simulation results. Indeed, the mea- 
surement errors in the independent variables can have signifi- 
cant impact on the regression analysis as also investigated by 
Kelly (2007), although the effect is only marginal in current 
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FIG. 2. — Comparison of five regression lines as afunction of the assigned 
amount of measurement errors on cr, using Graham et al. (201 1) sample. The 
percentages of assigned errors are given at lower right corners in each panel. 
In the case of measurement errors on independent variables above 10%, the 
difference between the fitted lines is clearly visible. 



datasets of the Mbh - c* relation. Typically the measurement 
errors in the velocity dispersion are assumed to be 5% or 10% 
in literature. 

In Figure|2]we compare the fitted lines to the Graham et al. 
Mbh~<^* dataset assuming measurement errors on <r* ranging 
from 0% to 30%. The difference between the fitted lines is no- 
ticeable and obvious when the measurement errors are large. 
Figure 3 compares the regression coefficients and intrinsic 
scatter derived from the five estimators as a function of the 
assigned errors. As a reference we show the results from the 
OLS estimator, i.e., for the unweighted fitting scheme without 
accounting for the intrinsic scatter as described in Isobe et al. 
(1990). This estimator is biased when there are measurement 
errors. 

For the intercept, the estimators do not give significantly 
different results except for the case of the BCES estimator. 
Both the intercept and slope estimated from the BCES estima- 
tor show very different behaviour in the high measurement er- 
ror regime, which is consistent with the result from Tremaine 
et al. (2002). 

For the slope, it is noted that all converge to the same 
value as the measurement errors in the velocity dispersion 
approach zero. Moreover, the estimators BCES, FITEXY, 
and Bayesian are very similar up to the 15% error level, 
thus indicating consistent estimation for these three estima- 
tors. The value of the slope from the BCES estimator be- 
comes higher compared to the others as the assumed errors 
on £r» increase. This is expected from the denominator term 
in Equation ((T). As can be seen, the estimated slope from the 
maximum likelihoodestimatoris almost identical to the 
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FIG. 3. — Direct comparison of regression results for the intercept (top), 
slope (middle), and intrinsic scatter (bottom) from the analysis of Figure ff] 
For the 5% errors the difference is only marginal. When the measurement er- 
rors are larger than 15% there is significant deviation between the estimators. 

OLS result up to measurement error amplitudes of ~ 16% on 
the independent variable. As noted and discussed in Kelly 
(2007) this biased behaviour is due to the implicit adoption of 
a naive uniform prior for the intrinsic distribution of the in- 
dependent variables; this bias is also noted in Kording et al 

(2006) . Based on this, we do not recommend the m aximum 
likelihood method as outlined in Section |2~3l It is sur- 
prising that the FITEXY with an ad hoc iterative approach 
gives fairly consistent results with that of the fully Bayesian 
approach (linmix_err). This is inconsistent with the result 
of the simulations performed by Kelly (2007). The source 
of this discrepancy is that for the FITEXY estimator Kelly 

(2007) did not refit the slope and intercept each time after the 
intrinsic scatter term is iteratively adjusted (see also, Kelly 
201 1). Instead, he just assigned the intrinsic scatter value such 
that the reduced \ 2 is equal to unity using the first minimiza- 
tion result of a and f3 for the zero intrinsic scatter case (i.e., 
just simply increasing a mt without re-optimization each time). 

For the intrinsic scatter, its level is very sensitive to 
the magnitude of the assumed measurement errors. Only 
1 i nmi x_e r r recovers a non-zero intrinsic scatter amplitude 
even in the case of assuming large measurement errors on cr*. 

3.3. Monte Carlo simulations 

An incorrect model for the distribution of the true values of 
x and y can lead to biased slope estimates, especially in the 
case of relatively large measurement errors on the indepen- 
dent variables, as pointed out by Gull (1989) and Kelly (2007) 
(see also, Auger et al. 2010, Mantz et al. 2010, and March et 
al. 2011). Here we use Monte Carlo simulations to inves- 
tigate the effect of an incorrect assumption on the intrinsic 
distribution of the independent variables. First, we generate 
three 10,000 simulated datasets by assuming respectively uni- 
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simulated sample: uniform distribution simulated sample: normal distribution 




7.9 8.0 8.1 4 5 6 0.0 0.2 0.4 0.6 7.9 8.0 8.1 4 5 6 0.0 0.2 0.4 0.6 

intercept slope intrinsic scatter intercept slope intrinsic scatter 



simulated sample: power-law distribution 




intercept slope intrinsic scatter 



FIG. 4. — Monte Carlo simulation results for the cases of the given uniform (upper left), normal (upper right), and power-law (bottom) distributions on x. 
Each column shows the distribution of intercept, slope, and intrinsic scatter estimated from the simulated datasets using the various estimators. FITEXYKdiy is 
the version of FITEXY estimator implemented by Kelly ( 2007 ). MLEjg means the maximum likelihood estimator with a single Gaussian model for the 
distribution of independent variable as described in Section |2~4l BPME3G is the Bayesian posterior median estimate using linmix_err procedure based on the 
normal mixture model with 3 Gaussians. The median value of the simulated distribution is plotted as a vertical solid black line while the true value is indicated 
by the red dashed line in each panel. 
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FIG. 5. — Histograms of mass residuals from the Mbh — o"* relation based 
on the FITEXY estimator. A single Gaussian fit with the center fixed to zero 
is expressed as a solid line. 

form, normal, and power-law distributions for x. The number 
of data points in each realization is set to be the same as that 
of the Graham Mbh - cr* data set (i.e., 64). The true intercept, 
slope, and standard deviation of Gaussian intrinsic scatter are 
assumed to be 8, 5, and 0.3 dex respectively, similar to typical 
values from the regression results given in Table Q] In other 
words, the sample points (y) from the given intrinsic relation 
(y = 8 + 5x) are scattered by the Gaussian random offsets with 
a = 0.3. Then Gaussian random noises, having zero mean and 
standard deviations equal to the measurement errors from the 
Graham et al. (201 1) data set, are added to both x and y. We 
fit the simulated data sets using the regression methods de- 
scribed in Section|2] 

Figure |4] shows the simulation results for the uniform, nor- 
mal, and power-law distributions of x, respectively. As al- 
ready pointed above, the estimated intercept and slope from 
MLE are biased and distributed similarly to that of the OLS 
estimator. This bias is regardless of the form for the intrin- 
sic distribution, and surprisingly the MLE still has a bias for 
the simulated sample from the uniform distribution. This is 
because the the maximum-likelihoodmethod assumes a 
uniform distribution for the independent variable in the range 
of —oo to oo, while in the simulations performed here the uni- 
formly distributed data have some finite range (i.e., fixed to 
be same as the range of Graham et al. data). As can be seen, 
the true values are well recovered if the likelihood function is 
changed to assuming a Gaussian for the distri bution of inde- 
pendent variables, as described in Section |2~4| (MLE i r, ) . This 
modified maximum likelihood estimates give very similar dis- 
tributions to that of the fully Bayesian estimates based on 
the linmix_err procedure (BPME3G). Here we also show 
the result of the version of FITEXY estimator used by Kelly 
(2007) to compare directly (FITEXYKeiiy). The biased be- 
haviour is same as noted in Kelly (2007). Thus this means that 
his implementation of FITEXY is inefficient compared to the 
one implemented here. From a viewpoint of how well the true 
values are recovered, all of BCES, FITEXY, and Bayesian 
estimators performed very well in this test. 

3.4. The sample difference 




6 7 8 9 10 

log (M m /M Q ) from Graham et al. 



O barred from McConnell et al. 
O barred from Graham et al. 




offset = -0.005 dex (-1%) 
scatter = 0.059 dex (15%) 



2.0 2.5 
log (<7,/km s~') from Graham et al. 

FIG. 6. — Comparison of the values of Mbh and a* for the galaxies in the 
set defined by the intersection of the McConnell et al. (201 1) and Graham et 
al. (2011) samples. The barred sample defined by McConnell et al. (2011) is 
marked with a red open square, while that defined by Graham et al. (2011) 
is marked with a blue open diamond. The dotted line indicates the identity 
relationship. 



In this section we investigate the sample discrepancy in de- 
tail using the most updated samples of McConnell et al. and 
Graham et al. We do this because the Mbh _ c* relations de- 
rived from these two datasets show a difference in the inter- 
cept (see Table Q]), which consequently affects the value of 
the virial factor. Note that the change in the M BH - c* rela- 
tion from the Giiltekin sample to the McConnell sample is 
obvious because ther e wa s a major update of the sample as 
discussed in Section 13.11 However, the difference between 
the McConnell and Graham data is not clear since these two 
sample have 50 galaxies in common. 

Figure [5] shows the mass residuals from the best-fit M B h - 
er* relation derived from the McConnell and Graham data. As 
can be seen, there is one extreme outlier, Circinus, in the Mc- 
Connell sample. Note that the central velocity dispersions of 
Circinus used in McConnell (i.e., 158 km/s) and Graham (i.e., 
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TABLE 2 

AGN BLACK HOLE MASSES AND STELLAR VELOCITY DISPERSIONS 



_i L 

2.0 



AGN sample 

_j L_ 

2.5 

log (ff,/km s -1 ) 



FIG . 7 . — The VPbh — cr* relation for the AGN sample (25). As shown in 
the right-hand side grey histogram, we are currently suffering from a lack of 
high mass AGN sample. Here we can clearly see the variability of the BCES 
estimator due to the effect of a single point (i.e., NGC 4253) that is subject to 
much larger measurement error than others. 

75 km/s) data are different from each other even though they 
were both taken from the HyperLEDAEI database. The cen- 
tral velocity dispersion for this object currently given in the 
HyperLEDA is 157.6 ± 18.8 km/s. We verified that the value 
listed by Graham et al. was a typo (Alister W. Graham 2012, 
private communication). 

In order to investigate the difference in the sample in com- 
mon between Graham and McConnell, in Figure [6] we com- 
pare the values of the black hole masses and stellar velocity 
dispersions. For the Mbh values, there are quite a few galaxies 
deviating from the identity relation. This is mostly due to re- 
cent updates of the black hole masses by Schulze & Gebhardt 
(2011) in the McConnell data. For the <r*, the values of the 
Graham sample are slightly larger on average than those of 
the McConnell sample, except for a few outliers. This slight 
average difference stems from the difference in the adopted 
velocity dispersion measures. Graham et al. (2011) used the 
central velocity dispersion provided in HyperLEDA, while 
McConnell et al. (2011) mostly used the effective velocity 
dispersion whenever it was available. This leads to systematic 
differences in the dispersion values as discussed in Tremaine 
et al. (2002). These mass and dispersion differences work 
to make the intercept smaller in the Graham sample com- 
pared to the McConnell sample. However, we note that the 
barred sample does not show any significant difference be- 
tween these datasets. We also performed the regression for 
the common sample only, and found that the intercepts re- 
mained almost the same while the slopes were reduced by 
0.2-0.3 compared to that of the entire sample. Therefore, the 
difference of the intercepts between McConnell and Graham 
samples is due to the different values adopted for the common 
sample, while the difference in slopes is primarily due to the 
non-overlapping sample. 

4. THE VIRIAL FACTOR 

The virial factor is of fundamental importance for estimat- 
ing AGN black hole masses in that it properly calibrates the 

5 http://leda.univ-lyonl.fr/ 
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NOTE. — Col. (1) name. Col. (2) virial product (VP B h =M m /f) based 
on the line dispersion (<7ij nc ) from reverberation mapping. Col. (3) reference 
for virial product. 1. Peterson et al. 2004; 2. Bentz et al. 2006b; 3. Denney et 
al. 2006; 4. Bentz et al. 2009b; 5. Denney et al. 2010; 6. Park et al. 2012; 7. 
Barth et al. 2011. Col. (4) stellar velocity dispersion. Col. (5) reference for 
stellar velocity dispersion. 1. Nelson et al. 2004; 2. Nelson & Whittle 1995; 
3. Ferrarese et al. 2001; 4. Onken et al. 2004; 5. Watson et al. 2008; 6. Woo 
etal. 2010; 7. Barth etal. 2011. 



measured virial product to a black hole mass for both the re- 
verberation mapping method and the single-epoch method. 
Following Onken et al. (2004), we determine the average 
virial factor (/) by forcing the AGN black hole masses onto 
the Mbh-ct, relation of quiescent galaxies. The AGN sample 
used here is listed in Table [2] with the corresponding refer- 
ences. We updated the AGN sample given in Table 2 of Woo 
et al. (2010) by updating the virial products from Denney 
et al. (2010), revising the rms line widths from Park et al. 
(2012), and including the new Mbh estimate for Mrk 50 from 
Barth et al. (2011). In Figure [7] we estimate the VPbh-c* 
relation with four regression methods, as in Figure Q] The re- 
gression results are listed in Table [3] The slope appears to be 
marginally lower than that for quiescent galaxies. This small 
difference in slopes might be due to noise and AGN selec- 
tion effects or it could be intrinsic, indicating a difference be- 
tween passive and active galaxies (see Greene & Ho 2006 and 
Schulze & Wisotzki 201 1). Note that the current sample is not 
representative of the AGN population since there is a deficit 
of high-mass AGNs, for which stellar velocity dispersion is 
extremely difficult to measure due to the overwhelming AGN 
luminosity. 

To determine the average virial factor we use the FITEXY 
estimator, fixing the intercept and slope to be the same as that 
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TABLE 3 

The VPbh-o"* Relation for the Active Galaxy 
Sample: log(VP BH /M ) = a + ,81og(<7*/200 km s" 1 ) 



Method 


a 


P 






Forward 


regression 




OLS 

MLE 

BCES 

FITEXY 

Bayesian 


7.25 ±0.14 
7.23 ±0.14 
7.30±0.17 

7.26 ±0.15 
7.24±0.17 


3.35±0.57 
3.20 ±0.59 
3.65 ±0.75 
3.38 ±0.63 
3.33 ±0.69 


0.41 ±0.06 
0.41 ±0.06 
0.43 ± 0.06 
0.47 ±0.09 




Inverse regression 




OLS 

MLE 

BCES 

FITEXY 

Bayesian 


7.74 ±0.23 
7.72 ±0.33 
7.70 ±0.30 
7.68 ±0.26 
7.68 ±0.37 


5.93 ±0.82 
5.88±1.21 
5.72±1.10 
5.68 ±0.95 
5.67 ± 1.87 


0.57 ±0.1 3 
0.51 ±0.12 
0.56 ±0.11 
0.62 ±0.28 



of the M B h - c* relation for the quiescent galaxies, 
2 _ (yi+iogf-a-pxj) 

X 2^ cr 2. + /3 2 cr 2. + 2 ■ W 
, =1 y,i ^ "x,/ Tu int 

Here, y = log(VP B H/M Q ) and x = log(cr»/200 km s" 1 ). The 
free parameters are only / and a mt . Adopting the regression 
results listed in Table Q] we estimate the virial factor and list 
the result in Table [4] 

Figure [8] shows the dependency of the virial factor on the 
adopted slope and intercept based on three datasets with four 
regression methods. As expected, the difference of / between 
the different regression methods for a particular dataset is 
small, with the only exception being the value of / obtained 
from MLE (red symbols). Estimated virial factors vary as 
much as a factor of 2 among the data sets, larger than the typi- 
cal uncertainties. The difference in / factors derived from the 
McConnell and Graham data is mostly due to the difference 
in the values from the sample of gala xies t hat overlap in these 
two data sets, as discussed in Section l3~4l The recent updates 
of Mbh measurements by Schulze & Gebhardt (201 1) lead to 
a smaller mean mass in the Graham sample compared to that 
in the McConnell sample. The difference of the adopted ve- 
locity dispersion measures results in a slightly larger velocity 
dispersion on average in the Graham sample compared to that 
of the McConnell sample. These combined differences make 
the intercept of the M B h _ cr* relationship smaller in the Gra- 
ham sample than in the McConnell sample, thus reducing the 
/ factor in the Graham sample regardless of the adopted re- 
gression methods. As can be seen, the derived virial factor 
is susceptible to the small variation of the quiescent galaxy 
Mftu-a* relation within the current calibration process. 

With the current AGN dataset, we are unable to constrain 
the / factor as a function of the mass range or host galaxy 
morphological type, since the number of sources in our sam- 
ple is small and the morphology of our sample is biased to- 
ward late-type galaxies. We note that a larger AGN sample 
(e.g., high-mass AGNs, especially) is needed for better statis- 
tical calibration. 

5. INVERSE FIT 

In addition to the conventional forward fit relation (i.e., 
fitting Mbh for a given c*, as we performed in previous 
sections), Graham et al. (2011) used an inverse fit for the 
Mbh-c* relation, suggesting that it corrects for possible sam- 
ple selection bias due to non-detection of intermediate-mass 
black holes (< 10 6 M Q ). We also follow their argument, re- 



TABLE4 

The derived average virial factor and intrinsic scatter 

BASED ON THE ADOPTED M BH ~<T* RELATION GIVEN IN TABLeQ] 



Method 

log{/) (7 mt 

From Giiltekin et 

licincr trvrwnrH rplntirvn 

MLE 0.82 ±0.09 0.43 ±0.05 
BCES 0.81 ±0.10 0.43±0.05 
FITEXY 0.81±0.10 0.43±0.05 
Bayesian 0.81 ±0.10 0.43±0.05 


l0g(/) <T mt 

al. (2009) Sample 

liQino" in\/prsp rpl ntioti 

0.55±0.12 0.54±0.06 
0.60 ±0.11 0.51 ±0.05 
0.59 ±0.11 0.52 ±0.06 
0.57±0.12 0.52±0.06 


From McConnell e 
using forward relation 

MLE 0.74 ±0.11 0.48 ±0.05 
BCES 0.72 ±0.11 0.49 ±0.05 
FITEXY 0.71 ±0.11 0.49 ±0.05 
Bayesian 0.71 ±0.11 0.49 ±0.05 


tal. (2011) Sample 

using inverse relation 

0.51 ±0.14 0.62±0.07 
0.52 ±0.13 0.62 ±0.07 
0.53 ±0.13 0.61 ±0.06 
0.52±0.13 0.61 ±0.07 


From Graham et 
using forward relation 

MLE 0.64 ±0.10 0.47 ±0.05 
BCES 0.55 ±0.11 0.50 ±0.05 
FITEXY 0.58 ±0.11 0.49 ±0.05 
Bayesian 0.58 ±0.11 0.49 ±0.05 


al. (2011) Sample 

using inverse relation 

0.42 ±0.1 3 0.58 ±0.06 
0.42 ±0.13 0.57 ±0.06 
0.45 ±0.12 0.56 ±0.06 
0.46 ±0.12 0.56 ±0.06 



fit all relations, and derive the average virial factors based on 
the inverse relations (see Table Q] [3] and HI). Note that ba- 
sically forward and inverse fittings are not the same in the 
presence of intrinsic scatter. Depending on the direction of 
regression (i.e., whether to choose Mbh or cr„ as the inde- 
pendent variable) the regressed slopes show large differences, 
leading to substantial changes in the virial factors. Therefore, 
we investigate and discuss the inverse fit in the context of the 
M B h _ relation. Now we have three factors related to the 
linear regression, which make the problem more complicated: 
measurements errors, intrinsic scatter, and truncation. 

If there is a truncation in the y-axis (i.e., logMBH) as argued 
by Graham et al. (2011), the conventional forward fit (fit y 
on x) causes a flattening in the estimated slope due to the in- 
creased loss of low mass black holes in the low cr* regime 
(e.g., see Appendix A in Mantz et al. 2010). The inverse fit 
(i.e., fit x on y) is not sensitive to this Malmquist-type bias, 
so long as incompleteness only exist in black hole mass. As 
shown by Kelly (2007), in order to avoid this selection bias 
on the regression result, it is necessary to assign the 'indepen- 
dent variable' as the variable used to select a sample. This 
approach has been generally adopted in the Tully-Fisher rela- 
tion studies since its sample is magnitude-selected and errors 
are smaller in magnitude than in velocity (e.g., Willick 1994; 
Tully & Pierce 2000; Bamford et al. 2006; Weiner et al. 2006; 
Koen & Lombard 2009; Williams et al. 2010; Miller et al. 
2011). 

In our sample, the measurement errors in the truncated co- 
ordinate (Mbh) are larger than in the other coordinate (<7«). 
Moreover, the sample selection might be highly inhomoge- 
neous and simple selection criteria may not be sufficient for 
describing it. The situation is more complex in the AGN 
sample selection. According to Schulze & Wisotzki (2011), 
even though the inverse relation is insensitive to the mass- 
dependent selection, it does not yield the intrinsic true rela- 
tion without incorporating the knowledge of the underlying 
host galaxy distribution function, which is currently hard to 
measure precisely. Furthermore, the AGN sample likely ex- 
hibits incompleteness in er* as well, as it is harder to measure 
cr* for AGN hosting more massive black holes due to their 
tendency to have higher luminosities. Thus there are good 
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FIG. 8. — Variation of the estimated virial factor as a function of the adopted intercept (left) and slope (right) of the quiescent galaxy Mbh relations 
taken from Table[T] Symbols mean the corresponding dataset used for estimation of the intercept and slope as expressed in upper right corner, while the colors 
of symbols indicate the regression methods used for them as given in the upper left corner. The horizontal dashed line indicates the value of the virial factor 
estimated from Woo et al. (2010) (i.e., 0.72). As an illustration, we add the green solid lines which show the dependence of the virial factor on the arbitrarily 
varied intercept (left) and slope (right) by fixing respectively the slope (left) and intercept (right) taken from the FITEXY estimates in the sample of McConnell 
etal. (2011). 



reasons to use either type of regression, but neither of them is 
completely free of bias. 

We provide both regression results in Table[T][3] and|4] In- 
verse regression results in a steeper slope compared to that 
of forward regression in the M BH - er* relation of quiescent 
galaxies. The calibration based on the inverse regression 
makes black hole masses inferred from the AGN virial prod- 
ucts smaller, since most of the AGN sample is located at the 
low-mass regime, thus leading to a reduction in the average 
virial factor. This biased dependency toward the low-mass 
regime motivates expansion in the dynamic range of sample 
of AGN containing both reverberation mapping data and mea- 
surement of er*. Based on these results, we conclude that the 
origin of the difference in the recently reported virial factors 
(Woo et al. 2010 based on forward regression vs. Graham et 
al. 201 1 based on inverse regression) is mostly due to the di- 
rection of regression adopted (i.e., whether Mbh is considered 
the independent or dependent variable), as well as the differ- 
ence in the samples used to calibrate the mass estimates. 

Feigelson & Babu (1992) suggested that we should choose 
the regression method for individual cases depending on the 
scientific question being investigated. Here the purpose of 
deriving the Mbh _ c* relation for local quiescent galaxies 
is to calibrate AGN black hole masses with determining the 
virial factor. By properly comparing the black hole masses 
of inactive galaxies to virial products of active galaxies, the 
average virial factor is constrained as discussed in the previ- 
ous section. Thus it is desirable to adopt the type of regres- 
sion which yields the relation that minimizes the scatter in the 
black hole mass estimates (Graham & Driver 2007). It is more 
common to adopt the host spheroidal quantity as the indepen- 
dent variable because the scaling relations are often used to 
infer black hole mass using the host spheroidal quantities as 
a proxy. Considering this, and the fact that the AGN sample 
likely suffers from Malmquist bias in both Mbh and cr*, we 
prefer the calibration from the traditional forward regression. 



6. DISCUSSION AND CONCLUSIONS 

We investigated the differences in the derived M B h - er* re- 
lation and virial factor using the recently compiled datasets 
of quiescent and active galaxies. The investigated possible 
origins of the difference are the fitting methodology and the 
sample difference. 

For the difference in regression methods, we utilized 
and compared four linear regression techniques: FITEXY, 
Bayesian, BCES, and Maximum likelihood. With the 
current level of measurement errors of the M B h - c* dataset, 
all estimators except for the maximum likelihood es- 
timator show good performance and consistent results with 
each other. There is no significant difference between the 
estimators. However, the assigned size of measurement er- 
rors on c* can have a significant impact on the regression re- 
sults, especially for the BCES and Maximum likelihood 
estimators. The Maximum likelihood method using an 
implicit assumption of a uniform distribution for the intrin- 
sic distribution of the independent variables introduces a bias 
which is clearly noticeable when the measurement errors on 
the independent variable are large (e.g., above 10% errors in 
the Graham sample as shown in Figure |2j. Without prop- 
erly accounting for the form of the intrinsic distribution of the 
independent variable, MLE estimates are very similar to the 
OLS results. Therefore we do not recommend this method for 
regression analysis in general. Of course for the Mbh _ cr» re- 
gression analysis the difference in the estimated slope is only 
marginal at the current adopted level of uncertainty on cr» 
(5%). The BCES estimator is also one of the good estimators 
based on the current measurement error level on cr*, although 
it may be problematic if the error is larger. Based on our sim- 
ulation results, the FITEXY estimator shows slightly better 
performance and the least-biased result compared to the other 
methods, although the others also perform well and the differ- 
ences are marginal. This is consistent with the result of Novak 
et al. (2006), although they did not provide an explicit imple- 
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mentation of all of the methods, nor a specific quantitative 
comparison. In general, we recommend both the FITEXY 
and Bayesian estimators, although the latter is computa- 
tionally more intensive, especially when the measurement er- 
rors are large. However, we note that the Bayesian estima- 
tor has the advantage over the method of FITEXY in that it 
calculates the full probability distribution function (i.e., pos- 
terior) of the parameters for the given data, and hence pro- 
vides well-defined and reliable parameter uncertainties. In 
addition, the Bayesian method can incorporate upper lim- 
its, as can the method of Giiltekin et al. (2009), whereas 
the FITEXY cannot. If we perform the regression using the 
Bayesian method, for the Giiltekin sample including up- 
per limits as well as secure measurements, the result changes 
from (a = 8. 19 ± 0.07, /3 = 4.04 ± 0.40, er int = 0.42 ± 0.05) 
to (a = 8.13 ± 0.07, (3 = 4.32 ± 0.38, cr int = 0.43 ± 0.05), 
thus log(/) correspondingly decreases from 0.81 ±0.10 to 
0.70 ±0.10. As discussed in Tremaine et al. (2002), accu- 
rate and consistent estimation of an individual stellar velocity 
dispersion with a correct measurement uncertainty is still re- 
quired and it will be an important factor for better constraining 
the Mbh - c„ relation and virial factor. 

The difference in sample is the most important factor con- 
tributing to the differences in derived Mbh relations. Giil- 
tekin et al. (2009) noted that the late-type galaxy and pseu- 
dobulge population in their sample is the source of the dif- 
ference in intrinsic scatter measurements by comparing their 
sample to that of Tremaine et al. (2002). Greene et al. (2010) 
found that the late-type low-mass galaxies show large scat- 
ter and are offset relative to the Mbh _ &* relation of elliptical 
galaxies using the sample of megamaser disk galaxies. By ex- 
tending the work of Graham et al. (2008), recently Graham 
et al. (2011) showed that the fraction of barred galaxies in 
their sample alters the Mbh _ c* relationship by dividing their 
sample into barred and non-barred galaxies. 

According to these previous studies, the Mbh _ c* relation 
seems to be not universal. It varies depending on the mass 
range and galaxy type. Correspondingly, the average / factor 
is also significantly affected by the sample population, since 
the intercept and slope from the quiescent galaxy Mbh re- 
lation are directly used in the calibration process. As investi- 
gated in this study, the differences in the adopted sample con- 
tribute to the change of the virial factor. Moreover, the direc- 
tion of regression (forward vs. inverse) causes further changes 
in the virial factor. We showed that the derived / factors vary 
as much as a factor of 2, which is from a combined effect of 
the sample and regression used. These differences could be 
thought of as an additional systematic uncertainty in the AGN 
black hole mass estimation via the current calibration process 
of the virial factor, since there is no obvious physical founda- 
tion for the selection of the appropriate sample and direction 
of regression. 

The true average / factor should not be changed by the host 
galaxy type since there should be no direct physical link be- 
tween the AGN BLR geometry and the global morphology of 
host galaxies. Unfortunately, the estimated average / factor 
may be subject to biases due to its calibration based on a sin- 
gle Mbh-c* relationship. However, since the current sample 
is not large enough to calibrate the virial factor as a function 
of galaxy type, it is better to use a single value of the mean 
/ factor for AGN Mbh estimation in order to avoid additional 
systematic errors until we get enough direct measurements of 
the structure of the BLR for an each individual AGN. We note 
that an alternative method to measure AGN black hole mass 
that derives the virial factor through BLR modeling has been 
recently developed and applied to the reverberation data (e.g., 
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FIG. 9. — The updated Mbh - c* relations for the inactive (black) and ac- 
tive (red) samples using the FITEXY estimator for the forward regression 
(upper) and inverse regression (lower). Shaded regions show the \u (68%) 
confidence intervals on the best-fit line. The inactive sample is from Mc- 
Connell et al. (201 1) and is the most recent one. The active sample is given 
in TablefS] 

Pancoast, Brewer, & Treu 2011; Brewer et al. 2011; Pan- 
coast et al. 2012). Given the uncertainties in the / factor, 
when investigating evolutionary trends in the Mbh _ rela- 
tion based on SE estimates, we recommend to use self con- 
sistent samples and techniques at different redshifts. In other 
words, one should measure the SE black hole masses con- 
sistently for AGN samples at different redshifts by using the 
cross-calibrated recipes based on the same virial factor. In 
this way the virial factor should be very similar for all sam- 
ples and cancel out in the determination of the evolution of 
logMfiH under the assumption that the virial factor is not a 
function of redshift (e.g., Woo et al. 2008). 

Finally, we present the updated Mbh _ cr* relation for local 
active galaxies based on the FITEXY estimator in Figure |9l 
where the forward (inverse) regression result is a = 7.97 ± 
0.14, P = 3.38±0.61, and a mt = 0.42±0.06 (a = 8.17 ±0.27, 
/3 = 5.47 ± 1.01, and a mt = 0.52 ±0.11). The AGN black 
hole masses were converted from the virial products using the 
virial factor log(f) = 0.71 ±0.11 (log</) = 0.53 ±0.13) de- 
rived in Section |4] From a methodological point of view, we 
prefer the forward regression as discussed in Section[5] Thus 
our preferred value for the virial factor is 0.71 based on the 
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preferred forward FITEXY/Bayesian estimation with the 
most recent sample (McConnell et al. 2011). This value is 
consistent with that of Woo et al. (2010) (i.e., 0.72) and dif- 
fers from that of Graham et al. (201 1) (i.e., 0.45) by - 0.26 
dex. The difference arises from the combination of sample 
differences and regression differences. It is worth noticing 
that the bottom panel of Figure [9] shows slightly better agree- 
ment between the non-AGN and AGN Mbh - c* relations, 
which may indicate that the inverse regression has less bias 
than the forward one and thus might be more reliable. How- 
ever, this conclusion only holds if we assume that the active 
and inactive galaxies share the same Mbh _ cr* relationship. 
Considering these issues, it is still not conclusive whether the 
inverse method is preferable with the current datasets owing 



to selection effects and limited dynamic range of the AGN 
sample. 
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